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The process of coherent creation of particle - hole excitations by an electric field in graphene is 
quantitatively described using a dynamic "first quantized" approach. We calculate the evolution of 
current density, number of pairs and energy in ballistic regime using the tight binding model. The 
series in electric field strength E up to third order in both DC and AC are calculated. We show how 
the physics far from the two Dirac points enters various physical quantities in linear response and 
how it is related to the chiral anomaly. The third harmonic generation and the imaginary part of 
conductivity are obtained. It is shown that at certain time scale tni tx E~^^^ the physical behaviour 
dramatically changes and the perturbation theory breaks down. Beyond the linear response physics 
is explored using an exact solution of the first quantized equations. While for small electric fields 
the I-V curve is linear characterized by the universal minimal resistivity a = n/2(e'^ /h), at t > t„i 
the conductivity grows fast. The copious pair creation (with rate E^^'^), analogous to Schwinger's 
electron - positron pair creation from vacuum in QED, leads to creation of the electron - hole plasma 
at ballistic times of order t„; . This process is terminated by a relaxational recombination. 

PACS numbers: 81.05.Uw 73.20.Mf 73.23.Ad 
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I. INTRODUCTION 



It has been demonstrated recently that a graphene 
sheet, especially one suspended on leads, is one of the 
purest electronic systems. It became increasingly evident 
that electronic mobility in graphene is extremely large ex- 
ceeding that in best semiconductor 2D systems The 
scattering of charge carriers in suspended graphene sam- 
ples of submicron length is so negligible that the trans- 
port is balhsticjl, The ballistic flight time can be 
estimated as ttai = L/vg , where Vg is the graphene veloc- 
ity characterizing the "relativistic" spectrum of graphene 
near Dirac points and L is the length of the sample. Dur- 
ing the ballistic flight conductivity calculated neglect- 
ing interactions with phonons, ripplons, disorder and be- 
tween electrons, etc., are consistent with experiments at 
least near the Dirac point at which no carriers are present 

i- 

In a simplified model of a single graphene sheet (ne- 
glecting scattering processes and electron interactions) 
the chemical potential is located right between the va- 
lence and conductance bands and the Fermi "surface" 
consists of two Dirac points of the Brillouin zone Q . A 
lot of effort has been devoted theoretically to the question 
of dc and ac transport in pure graphene due to the sur- 
prising fact that the conductivity is finite without any 
dissipation process present. A widely accepted "stan- 
dard" value of the "minimal dc conductivity" at zero 
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was calculated using a simplified Dirac model pMlOj by 
a variety of methods. Direct application of the the Kubo 
formula at zero frequency, disorder strength, temperature 
and chemical potential 0,0] utilizes certain "regulariza- 
tions" like the irj regulator which is interpreted as an 
" infinitesimal disorder" . The regulator is removed at the 
end of the calculation. A more customary application of 
the Kubo formula starts with finite frequency. As noted 
by Zicgler and others 0] the order of limits makes a dif- 
ference and several other values different from ai were 
provided for the same system. The standard value cri is 
obtained using a rather unorthodox procedure when the 
dc limit w — > is made before the zero disorder strength 
limit 77 — > is taken. If the order of limits is reversed, 
one obtains M: 
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When the limit is taken holding oj — -q one can even 
obtain a value of 0-3 = tt^ Q, thus solving the "miss- 
ing tt" problem (the same value was obtained very re- 
cently using yet another regularization in 11[). Indeed, 
at least early experiments on graphene sheets on Si sub- 
strates provided values roughly 3 times larger than ai 
12 1 . Recent experiments on suspended graphene Q 
demonstrated that the dc conductivity is lower, 1.7cti, 
as temperature is reduced to AK . Hence, while CT3 seems 
to be inconsistent with experiment, one still faces the 
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question of what is the proper theoretical value. Since 
the conductivity of clean graphenc in the infinite sample 
is a well defined physical quantity, there cannot be any 
ambiguity as to its theoretical value. Another theoretical 
approach to the problem was the use of Landauer formal- 
ism to the graphene sheet conductance. The conductivity 
appears as a limiting value of infinite width W in this ap- 
proach 0, and one obtains ci. Also resumming the 
series in disorder by a self consistent Born approxima- 
tion gives a I Qjwhile other resummations can lead to a 
different result |13j . 

In contrast to the controversy with respect to the dc 
conductivity, both the experimental and the theoretical 
situation for the ac conductivity in the high frequency 
limit is quite different. The theoretically predicted value 
in the Dirac model is (T2 independent of frequency under 
condition lj » T/h The Dirac model becomes 

inapplicable when lo is of order of 7//1 w 4 • \Q^^Hz or 
larger, where 7 is the hopping energy of graphene. It was 
shown theoretically using the tight binding model and ex- 
perimentally in (l5| that the optical conductivity at fre- 
quencies higher or of order 7/ft becomes slightly larger 
than 02- Moreover, in light transmittance measurements 
at frequencies down to 2.5 • IQ^^Hz it was found equal 
to (T2 within 4%. The tight binding model does not con- 
tain any other time scales capable of changing the lim- 
iting value of the ac conductivity all the way to — > 0. 
Therefore one would expect that the dc conductivity is 
a2 rather than cri. 

As is shown in the present paper, using the dynamical 
approach (a brief account of some results was published 
in [1, this is indeed the case. The basic physical 

effect of the electric field is a coherent creation of elec- 
tron - hole pairs, mostly, but not exclusively, near the 
two Dirac points. To effectively describe this process 
we develop a dynamical approach to charge transport in 
clean graphene using the "first quantized" approach to 
pair creation physics similar to that developed in rela- 
tivistic physics [17| to describe electron - positron pairs 
creation. To better visualize the phenomenon of resistiv- 
ity without dissipation directly at zero temperature, dop- 
ing, frequency and disorder we describe an experimental 
situation as closely as possible by calculating directly the 
time evolution of the electric current after switching on 
an electric field. In this way the use of a rather formal 
Kubo or Landauer formalism is avoided and as a result no 
regularizations are needed. Although we consider an infi- 
nite sample, the dynamical approach allows us to obtain 
qualitative results for finite samples by introducing time 
cutoffs like ballistic flight time. Various other factors 
determining transport can be conveniently characterized 
by time scales like the relaxation time for scattering of 
phonons or impurities. 

We show in detail that some aspects the linear response 
physics are not dominated by the two Dirac points of the 
Brillouin zone at which the spectrum is gaplcss. For ex- 



ample, large contributions (infinite, when the size of the 
Brillouin zone 27r/a (a = ZA is the lattice spacing), is 
being considered infinite) to the conductivity from the 
vicinity of the Dirac points are canceled by contributions 
from the region between them. This phenomenon is re- 
lated to the chiral anomaly in field theory [l^ . We anal- 
yse the use of massless Dirac (Weyl) model approxima- 
tion to the tight binding model and propose an effective 
chirally invariant regularization for it. 

In addition to the analysis of the linear response, we 
determin the range of applicability of the linear response 
approximation by both calculation of higher orders in 
(dc and ac) electric field and solving the model exactly 
in the nonlinear electromagnetic response regime. Only 
the zero chemical potential case (no net charge) is con- 
sidered. In this respect the work is complimentary to 
an earlier study by Mikhailov [l^ in which ballistic non- 
linear electromagnetic response to an ac field was calcu- 
lated at finite chemical potential using the Boltzmann 
equation within a semi-classical approach and major ef- 
fects we study were omitted. We first obtain pertur- 
bative corrections up to the third order in the electric 
field E. At this order qualitatively new phenomena like 
the third harmonic generation and the imaginary part of 
conductivity appear. The calculation of the corrections 
allows us to estimate the time scale at which perturba- 
tion theory breaks down. At this scale, tni oc i?^^/^, the 
physical behaviour is expected to change qualitatively. 
Therefore one has to resort to nonperturbative methods. 
Certain aspects of nonlinear ac response at zero chemical 
potential were also studied recently by Mishchenko [20l |. 
In his work disorder (taken into account phenomenologi- 
cally) dominates the purely ballistic effects by cutting off 
the ballistic times at relaxation time scale before is 
reached. 

Physics of the simplest tight binding model beyond 
the linear response is explored via an exact solution of 
the first quantized equations. It is a peculiar property 
of the tight binding model in the dynamical approach, 
that the exact solution of the equations for arbitrary 
momentum ky can be reduced to that for ky = 0; the 
constant electric field being in the y direction. More- 
over, the remaining equations, using Floquet theory, can 
then be reduced to a recursion relation. We use this 
property to effectively calculate the long ballistic flight 
evolution of various physical quantities like the current 
density and the number of created pairs. While for small 
ballistic times, t < tni, the conductivity settles at a2, 
at t > tni the current grows linearly. This increase can 
be explained using Schwinger's electron - positron pair 
creation mechanism p]| . The pair creation is a two di- 
mensional version of that in QED with the creation rate 
proportional to E^^-^ 2^. This, in absence of relaxation 



channels (for times below ttai = L/vg), leads to the cre- 
ation of a neutral electron - hole plasma at ballistic times 
of order tni. This process cannot continue forever and 
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is terminated by a relaxational recombination. The ap- 
plicability of Schwinger's formula for the electron - hole 
pairs creation rate was debated over a long time and we 
set the limitation on the applicability of this exact for- 
mula to graphene. 

The rest of the paper is organized as follows. The tight 
binding model and the dynamical approach are described 
in section II. The perturbation theory in electric field is 
described in section III. The minimal conductivity is ob- 
tained and the role of states beyond the neighbourhoods 
of the Dirac points (and their relation to " axial anomaly" 
and the Nielsen-Ninomya theorem) is clarified. The dy- 
namical linear response approach to the ac field is con- 
sidered. In Section IV perturbative corrections beyond 
linear response (up to third order in E) are calculated . 
The third harmonic and inductive contributions to con- 
ductivity are discussed. The exact solution for arbitrary 
field is constructed using the Floquet theory in Section 
V. It is used in Section VI to discuss the pair creation 
rate, conductivity and speculate about the electron - hole 
plasma. Finally Section VII contains discussion and con- 
clusions. 



II. THE MODEL AND THE DYNAMIC 
APPROACH TO ITS PHYSICAL PROPERTIES 



The dynamical "first quantized" approach. 

Since the crucial physical effect of the field is mainly 
a coherent creation of electron - hole pairs (though not 
always near the Dirac points, see below), a convenient 
formalism to describe the pair creation is the "first quan- 
tized" formulation described in detail in [l^. The spec- 
trum before the electric field is switched on is divided into 
positive and negative energy parts describing the valence 
and conduction band: 

iJk {E = 0) Wk = -£kUk; Wk = /e^) ' 



i?k [E ^0)vk = eki'k; 
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where ffk = |/ik|- A second quantized state is uniquely 
characterized by the first quantized amplitude. 
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which is a "spinor" in the sublattice space. It obeys the 
matrix Schroedinger equation in sublattice space 
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The tight binding model in an electric field 

Electrons in graphene are described sufficiently accu- 
rately for our purposes by the 2D tight binding model of 
nearest neighbour interactions Q . The Hamiltonian in k 
space is 



where p is defined in Eq.([5]). The initial condition corre- 
sponding to a second quantized state at T = in which 
all the negative energy states are occupied and all the 
positive energy states are empty is 
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H = 
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with 7 ~ 2.7eV being the hopping energy; Si = 
I (O, A/3)and 62,3 = f (±^,~^^ are the locations of 

nearest neighbours separated by distance a ~ 3A] and 
the pseudospin index A, B denotes two triangular sub- 
lattices. In most parts of the paper we keep the function 
/ik arbitrary. Let us first consider the system in a con- 
stant and homogeneous electric field along the y direction 
switched on at t = 0. It is described by the minimal sub- 
stitution, 



(5) 



with vector potential A = {0,—cEt) for t > (e > 0). 
Later it is generalized to more general fields including the 
ac field. 



A physical quantity is usually conveniently written in 
terms of ip- We will be interested mostly in the current 
density (multiplied by factor 2 due to spin) 



Jy (t) ^ -2e ^ 
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as well as the energy and the density of the electron - 
hole pairs. The summation is over the honeycomb-lattice 
Brillouin zone, see Fig. 1, in which two Brillouin zones 
are outlined. The energy of the electrons is changing due 
to the applied electric field (with no dissipation in the 
ballistic regime, see however Discussion) and is given by 
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The amplitude of lifting an electron into the conduction 
band (defined for the Hamiltonian without the electric 
field) is {ip it) |wk), where the scalar product is defined 
by {iplcj}) = ip'{(t>i + V'2'/'2j ^^'^ consequently the density 
of pairs reads 
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III. LINEAR RESPONSE, THE PSEUDO - 
DIFFUSIVE BEHAVIOUR AND THE PARITY 
ANOMALY 

Expansion in electric field 



Since the Hamiltonian depends on time via the vector 
potential that shifts the momentum p, a more useful defi- 
nition of the density of pairs taking into account the shift 
of the momentum will be given in Section III. 



Expanding in the dimensionless electric field strength 
(assumed for simplicity homogeneous and constant 
after switching on the field at i = 0), -^j. = 
giEk* ^yj^ _|_ Qj-^Q obtains for the first correction 

the following equation: 



Units and conventions 



The units of energy, time and distance are defined by 
the microscopic values 7, = (~ 2.5 • 10~^^s) and 
a, correspondingly. The unit of the electric field will be 
i?o = ^ — lO^^V/m, so that the dimensionless electric 
field is £" = E/Eq and the unit of conductivity will be 
e^/H. Effectively one can set 7 = fi, = a = e = 1. In 
these units the first quantized equation reads 



(18) 



Writing the correction as — ak^k + /Sk'f^k, it takes the 
form 



(19) 



idt'4>2 
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where the tight binding Hamiltonian matrix element hp 
takes the form 
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with initial conditions ai^ {t ^ 0) = 1 , (3-^ {t = 0) = 0. We 
use the abbreviations 



and 
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~£t for the 



where h — 2 cos (^) , andpy 
dc field. We will use extensively its expansion in powers 
of Ay-. 



The coefficient /3k denotes the amplitude of accumulation 
of electrons in the conduction band, whereas ak denotes 
the amplitude of remaining in the valence band. Solving 
the equations one obtains 




For example, the off - diagonal matrix element for the 
p component of the current Jy is —h!^, so that the current 
density, using the first quantized approach, is 

Jy--^i:^l(j!*^ '0^)V (17) 

The next two sections deal with the perturbative treat- 
ment of the electric field and later (after having shown 
that it fails at certain time scale) we will switch to a 
nonperturbative method. 



«k = -^iH; P^='-^{l-e-^^^^'-2^e^t). (22) 
2 4ek 

The expansion for the current, Ea. ([TT|) . in the electric 
field up to the first order contains the following momen- 
tum k contributions, Jy = X^kJk^ 

Jk-J^+J^+Jt (23) 
The zero order contribution is 

Jk = -4(4 'i^) = 2ek/k+ (24) 

where 

" 24 ■ 
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The correction gives the 'paramagnetic' and 'diamag- 
netic' contributions to the current densities: 



K^. 



(26) 



= 2£(ak-a^)£k/k =^(^k) [2£ki - sin (2eki)] 



(27) 



Both corrections for a specific momentum k diverge at 
large baUistic times as t, however one still has to integrate 
over the valence band momenta. 




Integration over momenta and physical 
interpretation of the "quasi -Ohmic" behaviour. 



To first order in electric field the current density is 



Jy — Jq -\- erf, 



(28) 



where the zero order contribution can be written as a 
derivative with respect to ky of a periodic function 



Jo 



k 



Jk 



(29) 



It vanishes upon integration, in accordance with the 
Bloch theorem, since the Brillouin zone can be chosen 
in such a way that it exhibits periodicity in the field (y) 
direction. For example, we can integrate over the rectan- 
gular area shown in Fig.l containing two Brillouin zones: 
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1/2 
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(30) 



-27r/3i/2 ^k"S ~ ^fe^,27r/3i/2 



It should be noted that / ^ 
^fc^ _27r/3i/2 = at every fc^., and due to the continuity 
of Ekeven at the Dirac points. Therefore one is left with 
the linear response. 

The conductivity can be divided into an apparently lin- 
early divergent part and a finite part, a (t) ~ (Jbz {t) + 
c^DP (t). As will be explained shortly, the contributions 
to the first term at large t stem mostly from the im- 
mediate neighbourhoods of the two Dirac points, while 
contributions to the second come from whole Brillouin 
zone. The part diverging linearly with time can be writ- 
ten (after some algebra) as a full derivative with respect 
to fc„ : 



(JBZ (t) = t 



h^h*^-hlh'A -2sl 



k 

-4t' 



k 



(31) 



FIG. 1: Fig.l Constant energy map of the lioneycomb-lattice 
Brillouin zone. The rectangle outlines two Brillouin zones, 
the area of integration in Eqs. pil and l31|l . The circles show 
the vicinity of the two Dirac points with radius A for the 
integration in Ea. (|34[) . 



This too vanishes upon integration over the Brillouin 
zone, since it is again an integral of a derivative of a 
periodic function, albeit this cancellation is nontrivial. 
In Fig. 2a the integrand of Eq. ipTj) is plotted within the 
integration area specified above, Ea. (P(I)) .The integrand 
is negative along the line connecting two Dirac points 
along the field direction, for example 



K. = 27r 



1 1 



KR = 27r - 



1 

71 



(32) 



(recall that a = 1 in our units), and positive elsewhere. 
In Fig. 2b several cross sections for various values 
are shown. One observes, that as k^ approaches 27r/3 
the integrand goes to — oo at ky = 27r/-\/3, corresponding 
to Dirac point K^. Consequently for all kx ^ 27r/3,the 
integrand is a regular function, and thus the integral over 
ky vanishes. At k^ = 27r/3 the integral over ky is finite, 
yet does not influence the two-dimensional integral. 

The remaining part 

aop {t) = -2 ^ (/k )' sin {2eut) , (33) 

k 

gives a finite result. Unlike the divergent part of the 
conductivity, the integral for t >> 1 (in units of tj) is 
dominated by contributions from the vicinity of the two 
Dirac points (circles in Fig. 2a of radius A ~ 1 (in units 
of h/a)). Indeed the prefactor {l^) is bounded from 
above by A^^ and integral over the area outside the cir- 
cles vanishes at i >> 1/A. The contribution of a single 
Dirac point is obtained by integrating to infinity (in polar 
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FIG. 2: Fig. 2a. The integrand of cjbz Eq.((3T| is plotted 
within the integration area shown in Fig. l.The integrand is 
negative along the line connecting two Dirac points along the 
field direction, and positive elsewhere. 




A physical picture of this resistivity without dissipa- 
tion is as follows. The electric field creates electron - 
hole excitations in the vicinity of the Dirac points in 
which electrons behave as massless relativistic fermions 
with the graphene velocity Vg playing the role of the light 
velocity. For such particles the absolute value of the ve- 
locity is Vg and cannot be altered by the electric field and 
is not related to the wave vector k. On the other hand, 
the orientation of the velocity is influenced by the applied 
field. The electric current is ev, thus depending on ori- 
entation, so that its projection on the field direction y is 
increased by the field. An important observation, made 
for example in ref.ji^lj is that the electric field, while 
creating charge transport, does not change the overall 
momentum. As a consequence the effects of impurities 
do not affect the pair creation in the same way they af- 
fect "free" carriers. These two sources of current, namely 
creation and reorientation are roughly of the same size 
in the linear response. As we will see below, the situa- 
tion changes at ballistic times when the linear response 
breaks down. 



Pair creation rate 

In analogy to the linear response in the current, the 
pair creation rate per unit area, as defined in Eq. ()13|) . 
can be calculated and to leading order in £ is: 



dt 



(36) 



FIG. 3; Fig. 2b. Several cross sections of the integrand of 
o'BZ Eq. (|31|I are shown for various values. 



coordinates centered at the Dirac point) 



C^DP 



fsinM^^^^il^ (34) 

(ZTT) Jip=oJq=0 Q 
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47r 



q'=0 



sin (q') 



q' 



At long time the upper limit can be replaced by infinity 
which results 



0'Z5P 



1 

2^ 



Si (oo) 



(35) 



Here one can safely multiply the result by the valley de- 
generacy 2. Returning to physical units, one obtains 
a = f72, Eq.@. 

The dependence on time was calculated numerically. 
After an initial increase on the natural time scale = 
h/j the conductivity approaches a2 via oscillations, see 
Fig. 1 in ref. [1] . The amplitude of oscillations decays as 



a power 



1 + 



i(2t) 



(for t » 



The result of the numerical integration over the Brillouin 
zone was given in Fig. 2 of ref. [l6| . It is dominated by 
the two Dirac points and at large times (compared to tj) 
behaves as 



at TT 



(37) 



It is however well known that a dc electric field in QED 
(even with massive relativistic electrons) renders the vac- 
uum unstable [171 121|. with the "rcnormalizcd" number 
of pairs depending on Therefore nonperturbative 

effects arc expected and will be studied below . 

The notion of " rcnormalizcd" number of pairs is a con- 
sequence of the fact that for unstable systems the Fermi 
level should be continuously " renormalized" as function 
of time. In the first quantized formalism this corre- 
sponds to a continuous modification of the " initial" state 
— > = Vk-ft defining holes. This leads to the fol- 
lowing definition of the pair density 37[ , 



NAt)=2 |(^(i)kp)|' = 2^ 



keBZ 
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(38) 
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It is used in the framework of the relativistic Dirac model 



and was recently extensively discussed in ref . [22| in con- 
nection with graphene (using the instanton approxima- 
tion). 



ac field 

A similar calculation within the dynamic first quanti- 
zation formalism for the evolution of the current density 
can be performed for an arbitrary time dependence of the 
homogeneous electric field. Let us consider an arbitrary 
time dependence of the y component of the vector poten- 
tial Ay = —£a (t) , subject to a " switching on" condition, 
a (t) = 0, for i < 0. Fourier transforms are defined by 



ait)= I *a (u) ; (0 - / e"^'^*? k W , (39) 



and similarly for the current density and other 
quantities. The next to leading order in £ first quantized 
tight binding equations are: 



Re a 



4 ^ {h*h' -h*'hf 

ill! ^— ^ p 



X lim 



s-)-0 



(44) 



Im a = 



4 ^ {h*h' - h*'hy 

if,} ^ — ^ p 



X lim 



2ujs 



s-j-O 



0. 



Therefore ac and dc conductivity are equal, as was also 
shown in recent calculations [1^. This is consistent with 
both the Kubo formula derivations and optical ex- 



periments jl5l | . Similar results were obtained in multiple 
layered graphene [24 1. 

A similar calculation can be performed directly in the t 
space (similar to the dc case) and shows that after a short 
transient one obtains the (T2 value for the ac conductivity 
independent of frequency. 



£k - W 

K 



Ek - 1^ 



ekH-0. 

(40) 

The switching on condition, ^-^ {t <{)) ~ 0, can be taken 
into account by the uj + iQ'^ substitution regularizing the 
way the Fourier transform is defined for zero frequency. 
From Ea.(|40| one obtains 



a {ijj) 



Ck H = - 



ujh*h' /e - h*h' - h*'h 
\/2a; (2£ - w) V ^^^*' + s^*' + h*'^h'/e 

(41) 

In the same order the ac conductivity is an integral over 
the sum of the "paramagnetic" and the " diamagnetic" 
contributions given by 



IV. BEYOND LINEAR RESPONSE. WEYL 
FERMIONS, PARITY ANOMALY AND THE 
NIELSEN - NINOMIYA THEOREM 

As will be discussed below, an electric field should not 
necessarily be very small in graphene experiments, so 
that corrections to linear response are of interest. In 
addition it is important to determine at what ballistic 
time scale perturbation theory breaks down. In this sec- 
tion we present a perturbative calculation of the leading 
nonlinear effect in both dc and ac, and discuss the ways 
to regularize the Weyl model "correctly" in order to cal- 
culate these corrections. 



Ultrarelativistic fermions near Dirac points 



jk(c^) {h*h'-h*'hY 



a{uj) £ (2e + w) (2e - w) 



h*h" + hh*" 
2 ^ . (42) 



Subtracting a full derivative (independent of frequency) , 
one can rewrite this as 



jk(cj) _ {h*h' -h*'hY u? 
a (w) £'^ 



4£2 - OJ^ 



(43) 



The tight binding model employed here has two Dirac 
points around which the spectrum becomes ultra - rel- 
ativistic e = Vg |AL__Rk|, Ai^ijk = k - see Fig.l. 



^. The ma- 



in our units the graphene velocity is Vg 
trix element of the tight binding Hamiltonian can be ex- 
panded around as 



^k = 



(Afc^ 4-iAfcj,). (45) 



Real and imaginary parts, taking the uj + iQ'^ definition 
into account, are 



The Weyl field describing the "left handed" [ij, l25| 
fermions ip^ , namely a field satisfying the relativistic 
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Dirac equation with zero fcrmion mass 
zc^tTAf = Vg (dx + idy) -02 

idt1p2 = Vg {da: " idy) "^f , 

can be constructed by the unitary transformation 



■01 = -01 , 02 



(46) 



(47) 



The matrix element around the second Dirac point Kj^ 
is different, 



= Vg {Akx — iAky) 



(48) 



and the Weyl fermions arc "right handed" particles that 
obey 



«at0f 
idtip2 



Vg {dx 
Vg {dx 



idy) ^2 

idy) tpf, 



(49) 



without any unitary transformation required. They be- 
long to a different representation of the 2+1 dimensional 
Lorentz group [26 1. 

This is not just a peculiarity of the model, but a 
rath er g eneral feature [6| of massless fermions on a lat- 
tice [ISl, [27 1 . It is well known that in order to get a 
massless spectrum of fermionic excitations with any ul- 
traviolet cutoff (hexagonal lattice is an example), they 
come in multiple locations on the Brillouin zone (species 
"doubling"). In the Hamiltonian formulation the mul- 
tiplicity is 2^, where D is the dimensionality of space, 
Z? = 2 in graphene. In addition, the graphene fermions 
are "staggered", meaning the spinor components "live" 
on different sublattices of the honeycomb lattice. This 
reduces the doubling to 2^~^ = 2. The doubling is in- 
timately linked to the parity (discrete chiral symmetry) 
[isl [26j . The two Dirac points have opposite chirality so 
that there is no " parity anomaly" . 

It is sometimes claimed in condensed matter literature 
that, at least while doing linear response, one can concen- 
trate on the two neighbourhoods of the Dirac points and 
neglect the rest of the Brillouin zone. Even more, often 
the calculation's result is just multiplied by the factor 2 
(the valley degeneracy) . Below we show that this is gen- 
erally not an accurate description of what happens. This 
is not an academic question, since the low energy Weyl 
theory is simpler than the tight binding model and will be 
used in the next section to extend the calculations into 
higher orders in the electric field. One therefore needs 
a proper ultraviolet regularization of the Weyl model. 
Similar regularization issues are well known in field the- 
ory whenever chiral anomaly is encountered. Roughly, 
a "correct" regularization should respect the chiral sym- 
metry leading to important constraints on the number 
and charges of massless fermions. Otherwise the model 
is called "anomalous" and results of pcrturbativc calcu- 
lations become arbitrary. The most famous example is 



the requirement that in each generation of elementary 
particles (Icptons and quarks) the sum of charges should 
be zero [28|. 

We therefore discuss in some detail the cancellation of 
infrared divergencies and the correct application of the 
"ultra - relativistic" approximation to the tight binding 
model. 



Cancellation of ultraviolet divergencies and the 
approximate chiral symmetry. 



The ultra - relativistic approximation, Eq. (|45p . fails 
when applied to the first term in the expression for con- 
ductivity asz given in Eq. pi[) . At first glance the in- 
tegral in Eg. pip is dominated by the two Dirac points 
since the integrand diverges there, sec Fig. 2. For the 
both (widely separated) regions of the Brillouin zone, Kl 
and Kfi, it has the same asymptotic form 



hh* -h*h' 
2^3 



hh* 



h*h" 



V3. 



lAkl' 



(50) 



The integral over the neighbourhood of each Dirac point 
converges in the "infrared" (here meaning k — )■ K^^^) 
due to the Jacobian |Ak|, but is linearly divergent in 
the "ultraviolet" and both have the same sign, sec Fig. 2. 
Hence the integration cannot be extended to infinity and 
the size of the Brillouin zone serves as a natural ultravi- 
olet cutoff. 

It is important to note that there is no cancellation of 
the divergence between the Dirac points since both have 
the same sign! The divergencies are however canceled 
by the contributions from regions of the Brillouin zone 
between the Dirac points in which the ultra - relativis- 
tic approximation is not valid. Therefore in the "ohmic" 
regime during ballistic times one is not allowed to ne- 
glect states far from the Dirac points and replacement 
by the Weyl equation is incorrect. Due to cancellation 
of the whole "divergent" term, Ea. ([3T|) . one can devise 
an appropriate regularization in the UV in which these 
contributions are cancelled and even generalize the pro- 
cedure to higher orders in the electric field. We pro- 
pose and use such a scheme below in section IVC. Simple 
recipes like the momentum cutoff regularization (circles 
in Fig.l) or giving the fermions an infinitesimal mass 
commonly used, may lead to unphysical results. As a 
consequence more sophisticated rcgularizations like the 
C- function regularization [llj, "dimensional regulariza- 
tion" [i^l and "Pauli - Villars" were developed for contin- 
uous field theories like the Weyl model. The tight binding 
model is very similar in this respect to the "lattice Hamil- 
tonian" regularizations of the field theory (Hamiltonian 
meant here with time kept continuous) and also satisfies 
the chiral invariancc criterion Il8i . 
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Nonlinear response in dc. Where does the linear 
response break down? 

Since the current density is an odd function of the elec- 
tric field, the leading nonlinear correction to it appears 
in the third order in the field. The calculation along the 
lines of subsection IIIA is quite straightforward albeit te- 
dious. One can use the Wcyl model instead of the tight 
binding model due to the following reason. It is well 
known in field theory that generally chiral anomaly ef- 
fects, including the ultraviolet divergencies discussed in 
section HID, appear only in one loop calculations [28j . 
We checked and found that indeed up to the third order 
the expression is finite in the ultraviolet. Within the dy- 
namical approach it is natural to perform the calculations 
without Fourier transforming into the uj space. 

Up to the third order the current density is 



1 



64 



(51) 



The correction therefore is growing as a large power of 
the ballistic time. It becomes as large as the leading order 
for t = 23/23-i/4£-i/2^^^ Hence the perturbation theory 
breaks down on the time scale of 



(52) 



As will be seen in Section V, this agrees well with the 
crossover time obtained from a nonperturbative calcula- 
tion. Of course this time should be larger than any other 
possible relaxation time (due to impurities, phonons etc.) 
present in the system. A similar expansion was obtained 
for the number of electron- hole pairs [3l ■ 

The result is the same for the tight binding and the 
corresponding Weyl model. The only issue would be the 
first order calculation within the Weyl scheme since it is 
divergent in the ultraviolet. In Wcyl model we used the 
momentum cutoff that will be described in section IVB. 



One observes that the perturbation theory is inapplicable 
for 



u'^ < £ and t > uj/£. 



(55) 



This is less restrictive compared to the dc condition, 
Eq.dia) for uj> £^/h-\ 

The present formulas can be compared with the results 
of the dynamical calculation beyond linear response by 
Mishchenko [l^ in which a phenomcnological model of 
relaxation was employed. For the inverse relaxation time 
r, his nonperturbative result is 



Red 



2vj£^ 



w2r2 



1/2 



1 - 



3f2 



16a;2r2 



0{£^ 



This is consistent with the second correction terms in 
Eq.dSlD. 

The inductive part was absent in linear response and 
now appears for the first time: 



Im cr 



t£' 



(56) 



Due to the two conditions, Eq. (|55p . it is much smaller 
than a 2- The third harmonic is generated with the real 
part 



1 



29a;4 



£' 



while the inductive part is absent at the present order. 



V. THE EXACT SOLUTION OF THE FIRST 
QUANTIZED TIGHT BINDING EQUATIONS 

Application of the Floquet theory and reduction to 

ID 



ac response and the third harmonic generation 

An analogous calculation for the ac field E = £ cos {cut) 
switched on at t = results in 



It is a peculiar property of the tight binding matrix 
Eq. (|4]) that the solution for arbitrary ky can be reduced to 
that for ky = 0. Shifting the time variable tot = t—ky/£, 
one can define an 'universal wave function' 



jit)/£ 



1 



-O {£^ 



£^ + Rca cos (uit) + Im cr sin {uit) + a^^ cos (Swt) 



^(i) =^(t) 



(57) 



T^i^^chroedinger equation (|14p is now void of any ky 
dependence. 



The first term is just the reflection of the initial condi- 
tions and is non - universal. In the limit w — > one 
recovers the dc result. The corrected value of the real 
(dissipative) part of the ac conductivity is 



-2int 



(58) 



Rea — - 
4 



63 
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V 128^4 54^2 



t'\£' 



(54) 



where the dimcnsionlcss frequency D, = £ / (2^/3) is in 
units of t~^. The ky component of the momentum enters 
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the solution via the initial conditions Eg. pl)]) only. In 
terms of the universal function it takes the form 



= 0) = {t = -ky/S) = uk. 



(59) 



Taking the Fourier transform of Eqs. ((58)) . the equa- 
tions in frequency space turn into recursion relations: 



(Uj) 
(w) 



--02 (uJ 
'"01 



2n) 
2n) 



H2 



n). 



(60) 
(61) 



Floquet theory (and the recursion relations) assure that 
the frequencies are discrete and form two series both in- 
dexed by an integer m 



(62) 



The "central" frequency v will take generally two in- 
commensurate values. Writing the Fourier amplitude 
ipi (wm) as pm, one obtains, combining the two equa- 
tions, Ea.(|5n| and Ea. (|5T|) . the forward and the backward 



recursions, 



UJrn — 



CJm — 5f2 (jJm — 217 



- 5n 



Prn~l 

(63) 



n 



in 



Pm+2, 



Pni+l 

(64) 



with normalization po ~ 1. The pi and have to be 
determined by the consistency between the forward and 
the backward equations. 

Generally there are two independent Floquet frequen- 
cies, however our system is special with the exact relation 



= 2f7, 



(65) 



which can be proven rigorously. This greatly simplifies 
the solution. To obtain a simple approximation, we ex- 
pand around the easily solvable case oi = n corre- 
sponding to 6 = 0. In this case the solution consists of 
two frequencies 



,±(0) 



(66) 



where the superscript denotes zeroth order in b. For 
experimentally accessible cases << 1, the Floquet fre- 
quencies are close to ±1. We expand the Fourier coeffi- 
cients in powers of 6: 



n(0) 



(67) 



where we do not distinguish for the time being between 
the two Floquet branches. Let us look for a solution with 
the initial choice 



P: 



(0) 



p 



(68) 



which in particular means that pi = at this order. It 
can be shown that the "central" frequency does not 
have a first order correction in b. Using this fact the 
forward recursion can be rewritten as: 

{ujjn + ^) (^m - 2fl0Jjn - 1) Pm (69) 
= b [(w„j + n)p,n-l + [uJrn - 2Jl)p,„+i] + b^ (w„i - 2VL) p^^. 

Inserting the expansions for pm and v into this relation 
yields the expression for v~^^'^\ 



M2) 



i/(o)/2-» 

^(0)2 - rj2 



/") - 50 



20) (21/ 



(0) 



m (2i/(o) -t- O) (1/(0)2 _ ^2^^ ■ 



50) 
(70) 



so that the positive Floquet frequency (up to second or- 
der) is = j^(o) _(_ 52^(2)^ Moreover, it turns out that 
the expansion in b converges in the whole relevant range, 
< 6 < 2, and greatly assists the numerical solution of 
the recursion relation. 

After the coefficients p^ and Floquet frequencies 
are found, the solution of Eq. ([58| is written in terms of 
the Fourier series: 



(71) 



s— ± m— — 00 

00 



s=± m= — 00 

The second line follows from Eq. (pT|) . ujf^ are defined in 
Egs. lp^ and (??) and we utilized Eq. lpS]) for simplifica- 
tion. The two complex coefficients are fixed by the 
initial conditions Eq. ([5^ . This solution is used to calcu- 
late the evolution of the current density, the energy and 
the number of electron - hole pairs. 



Results. Nonlinear regime and Bloch oscillations. 

The current density divided by the electric field (in 
physical units), a (t) = Jy (t) / E, is shown in Fig. 3 in 
ref. for various values of the dimensionless electric 
field E = E/Eq. The (microscopic) unit of the electric 
field is very large Eq = = lO^^V/m, so that at real- 
istic fields £ << 1. After an initial fast increase on the 
microscopic time scale tj, cr (t) approaches the universal 

2 

value o'2 = -It^ and settles there, consistent with the 
linear response, calculated in Section IIIB. Beyond the 
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0.5 - 
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12 3 4 

t in units of ty-\/ E / Eq 

FIG. 4: Fig. 3. The time evolution of tlie dc conductivity 
in units of (T2.The time is scaled with £^^^^ .The simulations 
were performed in the field range £ — 2^^^ — with UV 
cutoff equal to 1 /a. The straight line is the asymptotic linear 
behaviour given by EQ. (|72p . 










rn 


L=2*a 






1 \ L=2"a 


L=2''-''a 


M 




L=2« 


7 


// \ ^"^^^ 




1.0 












a 


if Landau. 


- Zener rate 


o 

0.8 




uniti 






c 






0.6 






■ol-S 







I 2 .3 4 5 



t in units of ty^J E / Eq 

FIG. 5: Fig. 4. The electron-hole pair creation rate as func- 
tion of time in units of £^^^^t^,ioi field £ = 2^^". The rate is 
scaled with £^^^. The ultraviolet cutoff is always l/a, while 
different curves represent different values for the infrared cut- 
off L/a. 



crossover time tni the conductivity rises linearly above 
the constant "universal" value a2'- 



J it) 
E 



(72 A£ 



1/2 



t 



(72) 



where A = 25/23i/*7r-2 ~ 0.75. This is consistent with 
the Landau - Zener calculation (instantons) within the 
Weyl model [s^, [sl], see below. The crossover time, de- 
fined by J [tni] /E ~ a2, is therefore 



7' 



(73) 



and is consistent with the perturbation theory breakdown 
ballistic time, Ea. ([5^ . This linear increase regime can be 
considered as a precursor of the Bloch oscillations, but 
is still universal in the sense that it depends solely on 
neighborhood of the Dirac points. 

The current on the time scale of order 



Att- 



eEva 



(74) 



exhibits Bloch oscillations, as is seen in Fig. 4 of ref. |16l |. 
where larger fields in the range £ = 2^* — are shown. 
It turns out that the current vanishes at points given ex- 
actly at multiples of ts/'i with being the period of 
the Bloch oscillations. One observes a peculiar feature 
that (apart from the " relativistic" initial constant seg- 
ment) the time dependence of a (t) is similar for different 
electric fields. Indeed, if one plots J/Ve versus tE, all 
the curves nearly coincide. Moreover, one obtains the 
following excellent fit for the current 



■m 

E 



-1/2. 



(75) 



The Bloch time is approximately the time required for 
the electric field to shift the momentum across the Bril- 
louin zone Apy = eEts ~ h/a. This time scale is very 
long for experimentally achieved fields, much longer than 
the ballistic fiight time. For example in a sample of sub- 
micron length, L = 0.5/Ltm, tbai — 2.3 • lO'^t^. If one 
assumes that the electric current can reach /max = 1mA, 
so for a typical width W = 1.5/im one has an elec- 
tric field -Einax = — lO^V/m corresponding to 
£ = 10^-^ (the voltage in such case would be quite large 



EraaxL ~ 5F). The first maximum of the Bloch 



oscillation will be seen at flight time of <b/4 = 3.6- 10 t^, 
which is of the same order as tf^ai ■ If one uses a value of 
the current typical to transport measurements / — 50/i^, 
the electric field is just E = 5-lO^V/m corresponding f = 
5-10-5 (voltage y = 250771^), ti3/4 = 7.2 • 10% » hai 
and is therefore out of reach. See, however a recent pro- 
posal 



Crossover at tni in the Weyl model 

It is expected that the transition to the nonlinear 
regime is dominated by the neighborhoods of the Dirac 
points. Therefore it can be obtained also within the Weyl 
model provided it is properly regularized in the ultravi- 
olet region consistent with chiral symmetry, as was dis- 
cussed in section IV. For the calculation of the current 
density it suffices to renormalize the electric current by 
subtracting the UV divergent term Eg. ((50)) : 



Jy = - 



4l'g 

{2ny 



|k|<A 



■rii^2) + ^t£ 



(76) 
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The results are practically indistinguishable compared 
to the tight binding model for times larger than a mi- 
croscopic time scale and much smaller than the Bloch 
time. Some examples of the tight binding model were 
presented in Fig. 3 of ref. In Fig. 3 the conduc- 

tivity in units of (T2 is given as a function of time scaled 
with E~^/^. The function is the same for all electric fields 
as can be seen from scaling properties of the Weyl equa- 
tions. The simulations were performed in the field range 
£ = 2-12-2"^ with UV cutoff equal to l/o. The straight 
line is the asymptotic linear behaviour given by Eq. (1721) . 
The asymptotic form is attained soon after tni , so it looks 
like quite a sharp crossover. 

Therefore this model can be effectively used to study 
the transition to the rapid creation of the electron - hole 
pairs and the creation of the electron - hole plasma, but 
naturally cannot be used to see the transition to the 
Bloch oscillation regime. 

Below we discuss the nature of the crossover to the non- 
linear behaviour and possible new physics of the electron 
- hole plasma which might emerge. 



VI. PHYSICAL PICTURE OF THE CROSSOVER 
TO THE NONLINEAR RESPONSE 

Where does the energy go? 

Beyond linear response one does not expect the cur- 
rent density to hold up to its linear response value in- 
definitely. The situation differs from that in diffusive 
systems in which energy dissipates into heat due to in- 
elastic scattering off impurities and phonons. In this case 
the system is not an isolated one but becomes part of a 
larger system including a " background" . In an isolated 
ballistic system in a constant electric field the energy ini- 
tially increases, as follows from the following argument. 
The total energy, Eq. ([T2|) . of electrons can be written in 
the first quantized formalism as 



Uit)^2{^it)\H\i;{t)). 



(77) 



It can be shown using Eq.Q, that the power in this 
driven system is proportional to the current density 



dt 



C/ = 2( V 



dt 



-H yj"^ = ^2eE (^ip ^ = EJy {t) , 

^ (78) 
like in dissipative systems. Consistently, as was shown 
above, the ballistic system has a finite conductivity like 
a dissipative system. We will refer to this as to "quasi - 
Ohmic" behaviour. 

Since the model does not provide a channel of dissipa- 
tion, where does this "Joule - like" heat crE^ go? The 
dynamical approach allows us to calculate the evolution 
of energy going beyond linear response. The energy of 



the system (calculated in a way similar to the current) 
is increasing continuously if no channel for dissipation is 
included. Therefore the conductivity originates in cre- 
ation of pairs near the Dirac points with an additional 
contribution due to the turning of particles toward the 
field's direction. To gain more insight into the nature of 
the crossover to nonlinear response we also calculated the 
evolution of number of the electron - hole pairs during the 
ballistic flight. 



Schwinger's pair creation formula and graphene 

The electron - hole creation rate by the dc electric field, 
-^Np, with rcnormalizcd Np defined in Ea. (|38|) . is shown 
in Fig. 4 as function of time for field £ = 2~^°. Such a 
field E = 2-^°Eo = lO'^V/m is quite realistic [U. The 
rate is scaled with £^/^, while the time is given in units 
of £~^^'^try. The ultraviolet cutoff is always 1/a, while 
different curves represent the infrared cutoff L/a. The 
length to width aspect ratio, L/W, is taken to be 1. The 
results do not depend significantly on it for all values 
presented, as long as 1/4 < L/W < 4. 

The time dependence of the rate exhibits the same time 
scales as the current density. At times smaller than t„/ 
the pcrturbativc formula, Eq. (j36[) is valid. Immediately 
after switching on the electric field (times of order t^) the 
rate behaves as t^. For < t < tni the pair creation rate 
per unit area rises linearly (with logarithmic corrections) 
according to Eq. ([57]) . However it is clear from Fig. 4 that 
the expansion breaks down at tni when the rate stabilizes, 
approaching the value (in our units of a~^t~^) 



1/2 



(79) 



The power dependence of the rate on electric field, 
i;3/2 is the same as the rate of the vacuum breakdown 
due to the pair production calculated beyond pertur- 
bation theory by Schwinger in the context of particle 
physics (when generalized to the 2-1-1 dimensions and 
zero fermion mass |2lJ, |22| ) . It is not surprising since the 
power 3/2 is dictated by the dimensionality, assuming 
the ultra - relativistic approximation is valid. However 
the physical meaning is somewhat different. 

Note that the definition of the (rcnormalizcd) parti- 
cle number, see Eq. (p8| . is different from the classical 
Schwinger's path integral definition (which actually de- 
termines the "vacuum decay" rate rater than the pair 
creation rate). The two are not the same within the Weyl 
model, as was shown asymptotically in the limit of large 
times in (STj . and we obtain their value. The calculation 
of the number of pairs can be approximately performed 
using the instanton approach initiated by Nussinov in the 
context of particle physics 37 1 and extended in (2^, [sl] . 
In condensed matter physics the method is known as the 
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Landau - Zcner probability j3Cl |36| . It provides an intu- 
itive picture of Schwinger's pair creation rate. Unfortu- 
nately this picture cannot be extended to ballistic times 
smaller than t^i in which one cannot use the asymptotic 
large time expressions. 

Note that in the Boltzmann equation approach to bal- 
listic transport developed in ref. 



231 the renormalized 



pair number was utilized. This is connected to a simple 
relation within the Weyl model between the rate and the 
current density, as was shown recently (30l |. 



VII. DISCUSSION, CONCLUSIONS AND 
GENERALIZATIONS 



To summarize, we studied the dynamics of the elec- 
tron - hole pair creation by an electric field in a single 
graphene sheet where the chemical potential is located 
right between the valence and the conduction bands. We 
assumed that the transport is purely ballistic, thus ne- 
glecting impurities, interaction with phonons, ripplons, 
as well as the screened interaction between electrons. The 
dynamical approach which originated in particle physics 
17| was adapted to the tight binding model of graphene 
and allowed us to calculate the evolution of the current 
density, energy and number of pairs beyond linear re- 
sponse. We clarified several delicate issues within lin- 
ear response like the correct dc conductivity value and 
a proper use of the ultra - relativistic Weyl model ap- 
proximation to the tight binding model. The question of 
proper regularization within the Weyl approximation was 
linked to the chiral (parity) symmetry and the anomaly 
cancellation in the tight binding model. Using proper 
regularization the leading correction to the linear re- 
sponse in both dc and ac fields was calculated. The ac 
response which is purely pseudo-dissipative (no inductive 
part) and frequency independent in linear response shows 
strong frequency dependence, third harmonic generation 
and the inductive behaviour. 

It should be emphasized that beyond linear response 
new scales appear. Generally the tight binding model 



in electric field has a time scale tj = h/^ and a di- 
mensionless parameter £ = E/Eq, Eq = 7/ (ea). In 
linear response the dimensionless parameter £ does not 
appear in conductivity. This explains why the ac con- 
ductivity is independent of frequency all the way from 
optical frequencies 1/t.y down to dc. Therefore it is not 
surprising that the dc conductivity is (T2 rather than ai, 
which was obtained in numerous calculations over the 
years d, 0, H, [l3|. Some papers, in which both values 
are obtained in different limits [8|, even raise the ques- 
tion, what quantity is actually measured in experiments 
in " bulk" graphene like those in ref. [l, Q ? Within the 
dynamical approach there is no room for ambiguity since 
no regularizations, limits or uncontrollable approxima- 
tions were used. Wc considered also finite size effects for 



periodic boundary conditions and found that the con- 
ductivity converges to (T2 for sizes of order W,L ^ 50a. 
We claim therefore that dc conductivity in graphene is 
a well defined physically measurable bulk quantity for a 
sufhciently large sample, W,L » a, and cannot have 
two different values within the same model. In particu- 
lar there is no dependence on the aspect ratio W/L for 
large samples contrary to the result of the Landauer ap- 
proach. Therefore it is important to ask what could go 
wrong in other approaches to the same quantity in the 
same model. 

The various approaches leading to this incorrect value 
(Ti can be broadly divided into two classes. One type of 
calculations involves as the first step the introduction of 
disorder as a way to regularize the problem. At a later 
stage the disorder strength is taken to zero 0, • The 
calculation generally involves an uncontrollable resum- 
mation of diagrams. In addition to the Lindhart dia- 
gram, other diagrams (infinite series) involving disorder 
are summed. For a "regular" system this leads to the 
Drude formula with the correct limit (infinite conductiv- 
ity) when the limit of vanishing disorder is taken. In 
graphene the nature of conductivity is different. There 
arc no charge carriers and the electric field first has to 
create the carriers (electron - hole pairs) and only then to 
accelerate them. The acceleration is also quite specific: 
the absolute value of velocity is fixed, while only the an- 
gle can change. The diagrams omitted in this process 
might not be small. There is no small parameter like the 
loffe - Regel h/epT, where r is the relaxation time, as in 
the "regular" case. One argues that due to "large" ep 
the other diagrams involving crossings are small. At the 
Dirac point ep = Q and this argument should be ques- 
tioned. Therefore the use of the simplest resummation 
might be the origin of the error. 

The second independent approach giving the same 
value (Ti originates in mcsoscopic physics. Disorder does 
not appear here and the only input is the description of 
"leads" and the boundary effects of the sample. Within 
the Landauer approach [10[ one counts the "evanescent" 
modes in a ribbon of finite width W and length L leading 
to a result (for the armchair boundary conditions) 

g2 

cr = 4— — Vcosh"2(7rnI^/i) ^ cti. (80) 

?i=0 ' 

It was claimed that The predictions including the finite 
sample conductivity and the Fano factor were confirmed 
by experiments [33| for very short samples on substrate. 
Note that the expression, Ea. ([5n)) . depends on the aspect 
ratio only, even when both W and L are large. We have 
performed calculations for finite samples with periodic 
boundary conditions (limited to 1/4 < Wj L < 4) and the 
results have a consistent large size limit of (T2 irrespective 
of the aspect ratio. Thus there is a discrepancy between 
the two methods which should be further clarified. 
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Beyond the hnear response the dimensionless param- 
eter £ in principle can give rise to "macroscopic" time 
scales t = S^tj. A rather unexpected result is that at a 
scale, 



eEvn 



(81) 



the physical behaviour qualitatively changes. This time 
scale becomes the same as the ballistic time thai = 



L/Vg 



2 • lQ^^t^ for length L = O.bfim for relatively 



weak fields £ = 10~^ corresponding to E = IQ^V/m. 
It is important to note that this is the only time scale 
that enters the Weyl model approximation to the tight 
binding. Within this model the microscopic scale does 
not appear and tni is the only combination of the avail- 
able parameters Vg and E. One could anticipate that the 
same scale appears in Weyl model as well [s^l. Larger 
scales however can be formed in the tight binding model. 
For example at a scale Ib = t^/£ Bloch oscillations set 
in. This is already beyond the Weyl approximation. 

We therefore summarize the behaviour of the conduc- 
tivity and the number of electron - hole pairs in the var- 
ious regimes at finite electric field E in turn. 

(i) Ballistic times t < tni- Pseudo - dissipative be- 
haviour. 

After a brief transient period (of order of several t^ = 
h/j) the current density of the tight binding model ap- 
proaches a finite value and stabilizes there. Therefore the 
minimal dc electric conductivity at zero temperature is 

2 

(72 = fx' '^^'^ pairs are created with various orienta- 
tions of velocity (the value of which is fixed at Vg) and 
part of the current (the polarization or " zitterbewegung" 
part) is due to reorientation of the charge carriers. In this 
regime the pair creation rate given by Eg. pTp is small 
enough so that pairs can be considered independent, pro- 
vided the number density of created pairs ^-p— is not 
too large. If it is very large the inverse process of pair 
annihilation (via various channels) cannot be neglected. 
Approaching t„i the perturbation theory breaks down. 
At this time the density of pairs is ^ • It becomes of the 
order of 10^^cm~^ for E = lO^V/m. Let us assume that 
this does not happen and proceed to the longer ballistic 
times. 

(i) Ballistic times tni < t < ts- Schwinger's pair cre- 
ation and formation of the electron - hole plasma. 

On the larger time scale nonperturbative methods 
should be used. The pair creation is more intensive and 
the pair density asymptotically follows an analog of the 
Schwinger's formula for massless fermions in 2-1-1 dimen- 
sions. Eg. ([M]) . The current above tni increases linearly 
with time 



J (t) = a2 



V3 



3/2 



E 



1/2 



In this regime most of the electrons are oriented along the 
field direction, so that J = 2eVgN (t). For example, in 
order to reach the density of = lO^^cm^^ in a sample 
with ballistic time tf, = L/vg = 2000t-y (L = 0.5/im) 
the field should reach E = lO'^V/m. If the sample is 
short enough the transport still can be ballistic, but due 
to its nonlinear nature a more likely scenario is that a 
dissipation channel opens up. 

One can only speculate what kind of dissipation pro- 
cess truncates the pair creation and stabilizes the electron 
- hole plasma created by the electric field. Of course, 
the standard candidates are collisions with impurities, 
phonons, ripplons and the electron - electron interactions. 
Here we point out that the system is "open" and one 
should consider the "radiation friction" scenario: pairs 
annihilate emitting photons which take energy out of the 
graphene sheet. The effects of radiation of energy into 
space might in principle be observable at elevated fields 
and should be investigated. 

Ballistic times t > ts- Bloch oscillations. 

If the system is still ballistic at yet longer times, Bloch 
oscillations set in with a period oits ~ ~t^/£. This 
time scale becomes the same as the ballistic time thai ~ 
2- 10^^^ for relatively weak fields £ = 10"'^ corresponding 
to E = IQ'^V/m, see Fig. 2 in ref. 0. While the 
Bloch oscillations are difficult to observe (see however a 
recent proposal (29|). the transition to a nonlinear regime 
is within reach of current experimental techniques. 

The dynamic approach was generalized to bilayer 
graphene (s^ l in which similar questions exist for a long 
time. The correct dc conductivity for the A^ layered 
graphene is equal to the dynamical one a ~ Na2, consis- 
tent with the vanishing frequency limit of the ac conduc- 
tivity The creation of the electron - hole plasma is 
even more likely in these systems. 

In this paper only ballistic transport was considered. 
In principle, disorder and electron - electron interactions 
could be incorporated within the dynamical app roach in 



app roac 



(82) 



the way the Boltzmann equation approach 
be extended beyond the linear response. In fact, using 
a phenomenological methodology, disorder has recently 
been incorporated for the ac field in ref. Similarly 
Coulomb interactions and the pair annihilations into pho- 
tons, phonons etc. can be taken into account. Generally 
though, when the system has a large number of elec- 
tron - hole pairs, the screening by the neutral plasma is 
more effective and the influence of impurities and inter- 
actions diminishes. Understanding these effects is crucial 
for investigating the (nonlinear) plasma waves and their 
damping 41 1. 

Let us note the relation between the dynamics on the 
time scale tni and quantum adiabatic transport near the 
quantum critical point 35, [s^. The calculation of the 
number of pairs can be approximately performed using 
the instanton approach initiated by in the context of par- 
ticle physics 37 1 and extended in (2^, [s^ . In condensed 
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matter physics the method is known as the Landau - 
Zcncr probabihty (sol . 36 1. It provides an intuitive pic- 
ture of the Schwinger's pair creation rate. The balhs- 
tic evolution in graphene therefore can be considered as 
an example of the adiabatic quantum evolution which 
attracted much attention recently in connection to the 
Kibble - Zurel mechanism of phase transition dynamics 
and others. Graphene dynamics at large ballistic times 
offers an accessible system in which these processes can 
be observed. 

Finally let us remark on the application of the dynam- 
ical approach to calculating the response to very short 
strong field pulses like the femtosecond (and an order of 
magnitude longer) laser pulses. A possibility of measur- 
ing the response to such fields was advocated by Rusin 
and Zawadzki ^3,- For this purpose the dynamical lin- 
ear response formulas, as presented in the Section III, can 
be directly applied for arbitrary time dependence of the 
pulse, while the nonlinear calculation of Section V de- 
scribes the step - like pulse of finite duration only. Other 
shapes of the pulse can be calculated by breaking the 
pulse into several segments of different constant field. 
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